

// includes, system
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <cuda.h>
#include <math.h>
// includes, project
#include <cutil_inline.h>

__global__ void bitSeqToCountArray(unsigned *a, unsigned int N, unsigned int *count_array, unsigned int nmer_size);
__global__ void andCountArrays(unsigned int *count_array1, unsigned int *count_array2,
			       unsigned int *count_array_and, unsigned int ca_length);
struct count_array {
	unsigned int nmer_size;
	unsigned int n_int_blocks;
	unsigned int * array;
};

/**
 * Copies the bitsequence b_h to the device b_d
 **/
 
extern "C" void bitsequenceToDevice(unsigned int *b_h, unsigned int *b_d, unsigned int size) {
	cudaMalloc((void **) &b_d, size);
	cudaMemcpy(b_d, b_h, size, cudaMemcpyHostToDevice);
}

void initCountArray_device(unsigned int nmer_size, unsigned int *count_array) {
	size_t size_ca = (pow(4,nmer_size))/8; //4^n bits, 8bits/byte
	cutilSafeCall(cudaMalloc((void **) &count_array, size_ca));
	cutilSafeCall(cudaMemset(count_array, 0, size_ca));
}

/**
 * Should take in an array of bitsequences
 * bitsequence to countarray (stores all countarrays)
 * delete bitsequence
 * @param nmer_size
 * @param bitsequence_array
 * @param numOfBitSeqs
 * @param bitsequence_sizes: array sizes of bitsequences
//TODO: compare arrays
//TODO: Figure correct grid/block/thread
//TODO: Implement parallel array reduction
 **/
extern "C" void compareSequences (unsigned int nmer_size, unsigned int ** bs_array,
				  unsigned int num_of_bs, unsigned int * bs_lengths) {
	
	unsigned int **ca_array_device;
	size_t size_ca = (pow(4,nmer_size))/8; //4^n bits, 8bits/byte
	int length_ca = size_ca/sizeof(int);
	// TEST: Print countarray length
	printf("%s%d\n", "Count Array Length: ", length_ca);
	// Array of CountArrays
	ca_array_device = (unsigned int **)malloc(num_of_bs*sizeof(unsigned int*));
	
	// Initialize distance matrix
	float dist_matrix[num_of_bs][num_of_bs];
	for (int i = 0; i < num_of_bs; i++) {
		for (int j = 0; j < num_of_bs; j++) {
			dist_matrix[i][j] = 0;
		}
	}
	// Convert all bitSeq to count_arrays =====================================
	for (int i = 0; i < num_of_bs; i++) {
		//for (int k = 0; k <10; k++) printf("%s%u\n","cu_bs: ", bs_array[i][k]);
		unsigned int numblocks = bs_lengths[i];
		
		// Allocate bitseq on device ===================================
		unsigned int *temp_bitseq; //Temp device bitsequence
		//printf("%s%u%s%u\n", "Size (numblocks) in Bitsequence", i, ": ", numblocks);
		printf("%s%u%s%u\n", "Size (numblocks) in Bitsequence", i, ": ", bs_lengths[i]);
		//size_t size = bs_sizes[i] * sizeof(int); //size of BitSeq at position i
		size_t size = bs_lengths[i] * sizeof(int);
		cutilSafeCall(cudaMalloc((void **) &temp_bitseq, size));
		cutilSafeCall(cudaMemcpy(temp_bitseq, bs_array[i], size, cudaMemcpyHostToDevice));
		
		// Allocate count_array on device ==============================
		cutilSafeCall(cudaMalloc((void **) &ca_array_device[i], size_ca));
		cutilSafeCall(cudaMemset(ca_array_device[i], 0, size_ca));
		
		// Calculate grid/block dimensions and run kernel===============
		int block_size = 384;
		int n_blocks = numblocks/block_size + (numblocks%block_size == 0 ? 0:1);
		bitSeqToCountArray <<< n_blocks, block_size >>> (temp_bitseq,
				numblocks, ca_array_device[i], nmer_size);
				
				
				
		/**TEST BITSEQ TO COUNTARRAY
		*unsigned int *ca_h1;
		*ca_h1 = (unsigned int *)malloc(size_ca);
		*cutilSafeCall(cudaMemcpy(ca_h1, ca_array_device[i], size_ca, cudaMemcpyDeviceToHost));
		*for (int q = 0; q < length_ca; q++) printf("%u",ca_h1[q]);
		**/		
		cutilSafeCall(cudaFree(temp_bitseq)); //TODO: necessary? cudaThreadExit?
	}
	
	//Compare all countarrays
	for (int i = 0; i < num_of_bs; i++) {
		for (int j = i+1; j < num_of_bs; j++) {
			unsigned int *result_ca_d;
			cutilSafeCall(cudaMalloc((void **) &result_ca_d, size_ca));
			cutilSafeCall(cudaMemset(result_ca_d, 0, size_ca));
			
			int block_size = 256;
			int n_blocks = length_ca/block_size + (length_ca%block_size == 0 ? 0:1);
			andCountArrays <<< n_blocks, block_size >>>
				(ca_array_device[i], ca_array_device[j], result_ca_d, length_ca);
				
			//Copy Countarray to host
			unsigned int *ca_h;
			ca_h = (unsigned int *)malloc(size_ca);
			cutilSafeCall(cudaMemcpy(ca_h, result_ca_d, size_ca, cudaMemcpyDeviceToHost));
			cutilSafeCall(cudaFree(result_ca_d));
			// Add up all 1 bits
			int bitcount = 0;
			for (int q = 0; q < length_ca; q++) {
				bitcount += ca_h[q];
			}
			//int distance = bitcount; //TODO: Calc actual distance
			float distance = 1.0 - (bitcount / min(bs_lengths[i], bs_lengths[j]))
			dist_matrix[i][j] = distance; 
			free(ca_h);
			
				
			/** TEST
			*copy countarray to host
			*unsigned int *ca_h;
			*ca_h = (unsigned int *)malloc(size_ca);
			*cutilSafeCall(cudaMemcpy(ca_h, result_ca_d, size_ca, cudaMemcpyDeviceToHost));
			*for (int p = 0; p < length_ca; p++) printf("%u",ca_h[p]);
			**/
				
		}
		
	}
	//Test print out dist matrix;
	for (int r = 0; r < num_of_bs; r++) {
		for (int s = 0; s < num_of_bs; s++) {
			printf("%f%s",dist_matrix[r][s], " ");
		}
		printf("\n");
	}
}

//TEST CODE
extern "C" void
runTest(unsigned int arraysize, unsigned int *b_h, unsigned int *b_d, unsigned int nmer_size)
{
	
	size_t size = arraysize * sizeof(int);
	
	
	cutilSafeCall(cudaMalloc((void **) &b_d, size));
	cutilSafeCall(cudaMemcpy(b_d, b_h, size, cudaMemcpyHostToDevice));
	
	//initialize a countarray
	size_t size_ca = (pow(4,nmer_size))/8; //4^n bits, 8bits/byte
	unsigned int * count_array;
	cutilSafeCall(cudaMalloc((void **) &count_array, size_ca));
	cutilSafeCall(cudaMemset(count_array, 0, size_ca));
	//initCountArray_device(nmer_size, count_array);
	
	
	
	int block_size = 128;
	int n_blocks = arraysize/block_size + (arraysize%block_size == 0 ? 0:1);
	bitSeqToCountArray <<< n_blocks, block_size >>> (b_d, arraysize, count_array, nmer_size);

	//copy countarray to host
	unsigned int *ca_h;
	ca_h = (unsigned int *)malloc(size_ca);
	cutilSafeCall(cudaMemcpy(ca_h, count_array, size_ca, cudaMemcpyDeviceToHost));
	for (int j = 0; j < (size_ca/4); j++) printf("%d%s%u\n", j,"CA: ",ca_h[j]);


	cutilSafeCall(cudaMemcpy(b_h, b_d, sizeof(int)*arraysize, cudaMemcpyDeviceToHost));
  	// Print results
  	for (int i=0; i<11; i++) printf("%u %u\n", i, b_h[i]);
  	// Cleanup
  	//free(b_h); 
	cudaFree(b_d);
 
    cudaThreadExit();
}







//=======================================================================================
//                      DEVICE CODE 
//=======================================================================================
// Chops bits to first n elements (for now)
// TODO: Sliding window across intblocks
__device__ unsigned int chopToNmer(unsigned int sequenceBlock, unsigned int nmer_size,
				   unsigned int window_pos ) {
	int shift_l = window_pos*2; //each element is 2 bits
	int shift_r = (sizeof(unsigned int) * 8) - (nmer_size * 2);
	sequenceBlock <<= shift_l;
	return sequenceBlock >>= shift_r;
}

// Takes in an nmer value and sets a 1 bit on the correct integer in the array
//TODO: Figure out __fdividef()
//TODO: BUG: some Countarray ints set to all 1's (4294967295)
__device__ void setBit(unsigned int nmer, unsigned int * count_array){
		unsigned int int_idx = nmer/(sizeof(int)*8);
		unsigned int bit_idx = nmer % (sizeof(int)*8);
		// Translates flag_pos to decimal ( pos 4 -> 0b1000 -> 8 )
		//unsigned int flag_pos_dec = __powf(2, (bit_idx - 1));
		unsigned int flag_pos_dec = 1 << (bit_idx -1);
		count_array[int_idx] |= flag_pos_dec;
	//}
}

/**
 * Converts bitsequence array to countarray by setting bit flags in
 * correct positions within an integer bit representation
 **/
__global__ void bitSeqToCountArray(unsigned int *a, unsigned int arraysize, 
				   unsigned int *count_array, unsigned int nmer_size)
{
	int idx = blockIdx.x * blockDim.x + threadIdx.x;
	int max_reads_per_block = (sizeof(int)*8) - nmer_size + 1; //TODO: move this out of kernel
	
	if (idx < arraysize) {
		for (int i = 0; i < max_reads_per_block; i++) {
			int nmer = chopToNmer(a[idx],nmer_size, i);
			setBit(nmer, count_array);
		}
		
		
	}
	__syncthreads(); 
}

/**
 * Multithreaded AND operation on two arrays
 * Calculates the Hamming weight (pop count) of each integer and
 * assigns the value back into the ANDed countarray
 **/
__global__ void andCountArrays(unsigned int *count_array1, unsigned int *count_array2,
			       unsigned int *count_array_and, unsigned int ca_length) {
	int idx = blockIdx.x * blockDim.x + threadIdx.x;
	if (idx < ca_length) {
		count_array_and[idx] = count_array1[idx] & count_array2[idx];
	}
	//Count bits in each integer in ANDed array
	int count = 0;
	if (idx < ca_length) {
		while (count_array_and[idx]!=0) {
			//if(count_array_and[idx] & 
			count += count_array_and[idx] &= 1;
			count_array_and[idx] >>= 1;
		}
		count_array_and[idx] = count;
	}
}




